Atomic collision dynamics in optical lattices 
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We simulate collisions between two atoms, which move in an optical lattice under the dipole- 
dipole interaction. The model describes simultaneously the two basic dynamical processes, namely 
the Sisyphus cooling of single atoms, and the light-induced inelastic collisions between them. We 
consider the J = 1/2 — » J = 3/2 laser cooling transition for Cs, Rb and Na. We find that the hotter 
atoms in a thermal sample are selectively lost or heated by the collisions, which modifies the steady 
state distribution of atomic velocities, reminiscent of the evaporative cooling process. 
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I. INTRODUCTION 

Laser cooling and trapping techniques have made it 
possible to study and manipulate samples of cold neutral 
atoms ffi] . It has lead to the precision control of atomic 
matter, and also opened new possibilities for investiga- 
tions of interaction dynamics between atoms, especially 
those mediated by light g. By controlling the inter- 
nal states of the atoms we obtain access to their center- 
of-mass motion |l|J^]. Sisyphus cooling and polarization 
gradient cooling methods 0, which allow one to break 
the Doppler limit for low temperatures, are based on the 
way polarized light connects to the atomic states. By 
combining more than one laser beam we can introduce a 
spatially changing polarization state into the total light 
field interacting with atoms. As light introduces Stark 
shifts on atomic states, the spatially changing polariza- 
tion appears as a spatially changing potential for atoms 
with a suitable angular momentum state structure. In 
addition to the cooling effect, this has made it possible 
to build periodic or quasiperiodic lattices |pj, where the 
light field traps the atoms. We describe the Sisyphus 
cooling and lattice structure in Sec. H. 

The basic laser cooling method is Doppler cooling JD, 
which is produced by the random scattering of photons 
absorbed from the laser beam. This cooling mechanism 
is not related to the polarization states of light. For 
alkali atoms this method has a limiting temperature, 
the Doppler temperature To — hT/2kB, where T is the 
atomic linewidth. Using the polarization states, i.e., Sisy- 
phus cooling and polarization gradient cooling, one can 
go below Try until the photon recoil limit, Tr, is reached, 
and creating a lattice as a byproduct. The values of Tr 
and To for used elements are given in Table |. Typically 
one reaches a thermal equilibrium where the atoms are 
more or less localized at lattice sites, but can also move 
between them. The efficiency and the degree of localiza- 
tion have been studied thoroughly in the past ||. Since 
the best filling ratios (number of atoms per site) with 
small-detuning lattices are on the order of 10 % H, one 
can consider the gas sample as noninteracting. Larger 



filling ratios have been achieved lately by special tech- 
niques in far-detuned optical lattices (6). Based on the 
experience in standard magneto-optical traps (MOTs), 
the increasing density for small-detuning lattices is ex- 
pected to lead to strong loss and heating of atoms due 
to collisions, which become strongly inelastic in the pres- 
ence of the near-resonant cooling light [g,f7| . By a "small- 
detuning lattice" we mean one that is detuned a few 
atomic linewidths below the transition frequency. 

The light-assisted collisions are based on the fact that 
the two slowly approaching atoms form a quasimolecule, 
which the light can excite resonantly during the ap- 
proach, even though the cooling beams are clearly off 
the resonance with the atomic transition energy which 
corresponds only asymptotically to the molecular state 
transition energy. The resonance occurs at relatively 
long distances jq|, where the dominant contribution to 
the molecular behavior, i.e., the interaction potential be- 
tween atoms, comes from the dipole-dipol e in teraction 
(DDI). We describe these collisions in Sec. 



rive an expression for DDI in Sec. IV. We have presented 



III, and de- 



the first results of our work in a previous short publica- 
tion ||, where details of the derivation were omitted, so 
here we present them in full. It should be noted that one 
can develop e.g. a mean field approach to atom-atom pro- 
cesses via the DDI [ITy 13 . In our approach, albeit with 
some limitations, we allow the atoms to move. Further- 
more, we consider the problem in the atom-atom basis, 
with the full Zeeman substate structure, in the presence 
of the cooling/lattice-building laser beams with spatially 
changing polarization structure. Thus our model treats 
the Sisyphus cooling and localization of the atoms, and 
the atomic collisions dynamics consistently, within the 
same framework. 

In order to describe two multistate atoms moving 
quantum mechanically as wave packets in position and 
momentum space, and being coupled both to the spa- 
tially changing laser field as well as the vacuum field 
producing spontaneous emission, we would in principle 
have to use the density matrix description. This is not 
computationally possible currently, but we can go around 



the problem by describing the spontaneous emission with 
quantum jumps. As described in Sec. fV], we use the 
Monte Carlo wave function (MCWF) method to build a 
statistical ensemble of time evolution histories, which ap- 
proximates adequately the actual density matrix pl|-fi6[ . 
To implement this approach on our description of atomic 
collisions in lattices is not straightforward, and in Sec. VI 
we describe the details related to numerical simulations. 
Inelastic collisions in MOTs have been modelled exten- 
sively with semiclassical models [JL7| . We describe these 
models briefly in Sec. VII. They provide a tool for under- 



standing some of the physics behind the numerical data, 
and for estimating the processes affecting the multitude 
of boundaries for the numerical methods. 

O ur M onte Carlo simulation results, presented in 
Sec. VIII, indicate that the hotter atoms in our ther- 



mal sample, due to their stronger mobility between the 
lattice sites, are more likely to collide inelastically than 
the colder ones. On the other hand, the simulations, sup- 
ported by semiclassical estimates, show that inelasticity 
plays a relevant role only if the atoms end up in the same 
lattice site simultaneously. In other words, the effect of 
the dipole-dipole interaction remains small if the atoms 
do not share the same lattice site. This, of course, also 
depends on the chosen laser field parameters such as in- 
tensity and detuning. When the close, same-site encouter 
occurs, however, it is most likely a strongly inelastic one 
leading to the loss of atoms. Basically, we see a process 
similar to evaporative cooling, where the hotter atoms are 
selectively heated or ejected from the lattice. It should be 
noted, though, that despite the fully quantum mechanical 
nature of our approach, it still omits many other effects 
affecting the atomic cloud in the lattice. Photons scat- 
tered incoherently by atoms can be reabsorbed, which 
produces a radiation pressure; this process also heats the 
atomic cloud as its density and thus optical thickness 
increases (j. 

The observations made in this paper follow those from 
our previous study |9[. The results given in this paper, 
however, have been obtained with an improved approach 
compared to Ref . |9[| , and we have also extended our stud- 
ies to all the basic alkali species. Furthermore, here we 
give the detailed description of our approach and its com- 
putational aspects. Finally, the discussion in Sec. IX con- 
cludes our presentation. 



II. SISYPHUS COOLING AND OPTICAL 
LATTICES 

In this Section we present the atom-laser system under 
study and describe briefly the basics of Sisyphus laser 
cooling of neutral atoms in an optical lattice. Detailed 
review of the subject can be found in Refs. 
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FIG. 1. The level structure of a single atom. We show the 
squares of the Clebsch-Gordan coefficients of corresponding 
transitions describing the strengths of couplings between the 
Zeeman sublevels. The difference between the laser frequency 
oil and the atomic frequency ujq, i.e., the laser detuning, is 5. 



A. Atom— laser system 

We consider here atoms having ground state angular 
momentum J g = 1/2 and excited state angular momen- 
tum J e = 3/2 corresponding to alkali metal elements 
when the hyperfme structure is neglected. The resonance 
frequency between the states is loq so that Tlluq = E e — E g 
where E e and E g are energies of the ground and the ex- 
cited states in zero field. A single atom has two ground 
state sublevels \g±\/2 > and four excited state sublevels 
l e ±3/2 > an d |e±i/2 > where the half-integer subscripts 
indicate the quantum number m of the angular momen- 
tum along the z direction, see Fig. 111. The values of 
atomic masses that are used in our simulations are those 
of cesium ( 133 Cs), rubidium ( 85 Rb) and sodium ( 23 Na) 
which are the conventional alkali metal elements used for 
laser cooling of neutral atoms, see Table |. 

The laser field consists of two counter-propagating 
beams with orthogonal linear polarizations and with fre- 
quency ujl- The total field has a polarization gradient in 
one dimension and reads 



E(z, t) =£ (e x e 



ik r z 



ie v e- lk >- z )e~ luJLt + c.c. 



(1) 



where £q is the amplitude and k r the wavenumber. With 
this field, the polarization changes from circular a~ to 
linear and back to circular in the opposite direction a + 
when z changes by Al/4 where Al is the wavelength of 
the lasers. 

The periodic polarization gradient of the laser field is 
reflected in the periodic light shifts, AC-Stark shifts, of 
the atomic sublevels creating the optical lattice structure. 
The relative strengths of the couplings between a single 
ground state sublevel and various excited state sublevels 
vary spatially according to the polarization of the light 
field due to unequal values of the Clebsch-Gordan coef- 
ficients for different transitions. This induces light shifts 
and produces a periodic optical potential structure such 
that the shape of the light-induced potentials is the same 
for the two ground state sublevels but the potentials are 
shifted spatially with respect to each other by Al/4, see 
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FIG. 2. Schematic view of the optical potentials for the 
two ground state Zeeman sublevels. Lattice structure is cre- 
ated due to the periodic polarization gradient of the laser 
field. 



Fig. 0. The top of the optical potential for one sublevel 
coincide with the bottom of the other one. 



B. Sisyphus cooling 

When the atomic motion occurs in a suitable velocity 
range, optical pumping of the atom between ground state 
sublevels reduces the kinetic energy of the atom Q . This 
occurs because within the suitable velocity range, quan- 
tum jumps and optical pumping to another ground state 
sublevel tend to occur when the atom is near the top of 
the optical potential and is transferred to the bottom of 
the other one due to a quantum jump. Thus the sub- 
sequently emitted photon has a larger energy than the 
absorbed one and the kinetic energy of the atom is there- 
fore reduced, and the atom is cooled. After several such 
cooling cycles the atom localizes into the optical poten- 
tial well, i.e., into an optical lattice site. Figured shows 
the optical pumping cycles between the ground state sub- 
levels cooling an atom, and the oscillations of the atomic 
wave packet after localization into an optical lattice site. 

The intensity of the laser field and the strength of the 
coupling between the field and the atom is described by 
the Rabi frequency f2 = 2d£$/Ti where d is the atomic 
dipole moment of the strongest transition between the 
ground and excited states. The detuning of the laser 
field from the atomic resonance is given by 5 = lul — ^o- 
As a unit for fi and S we use the atomic linewidth T. 



C. Localization in lattice 

When the steady state is reached after a certain period 
of cooling, atoms are to a large extent localized into the 
lattice sites. In this study we deal with near-resonant 
bright optical lattices where the laser field is detuned a 
few atomic linewidths to the red of the atomic transition. 
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FIG. 3. Sisyphus cooling and the localization of an atom 
into the optical lattice. We show a possible time evolution for 
a single atom wave packet for two ground state Zeeman levels, 
(a) m g = —1/2, (b) m g = +1/2. The result shows the optical 
pumping cycles and the localization of a single atom into the 
optical lattice. This example forms one member of a Monte 
Carlo simulation, and the discontinuous changes between the 
two ground states are due to quantum jump events from the 
excited state (not shown), selected to happen randomly with 
an appropriately weighted probability. If the run is repeated, 
the jumps would appear at different times again. 

TABLE I. Atomic properties. Masses M in a.u. and 
linewidth energies e — {%T)/E r (for the definition of the recoil 
unit E T see Table |l|) . The Doppler temperature Tn and the 
recoil temperature Tr = (h 2 k^\ /MIcb in fj,K. Here k r is the 
wavenumber of the laser. 



property 


Cs 


Rb 


Na 


M 


133 


85 


23 


£ 


2400 


1600 


400 


T D 


120 


142 


238 


Tr 


0.20 


0.37 


2.4 



The laser parameters 51 and 8 determine if the lattice is 
in the "jumping" or in the "oscillating" regime, depend- 
ing on the average number of atomic oscillations in a 
single lattice site before the atom is optically pumped to 
neighboring sites jig . It must be noted that tight local- 
ization and occupation of the lowest vibrational levels of 
a periodic lattice potential increases the optical pumping 
time t p and the time of localization within a single lattice 
site becomes longer compared to the semiclassical values 
presented in Table y. 

We are interested in the effect of inelastic collisions be- 
tween atoms in the presence of near-resonant light H in 
optical lattices. These collisions occur when two atoms 
occupy the same lattice site. To observe efficiently the ef- 
fect of inelastic collisions we have chosen $1 and S in most 
of the simulations so that the lattice is in the "jumping" 
regime, i.e., semiclassically speaking the atoms on aver- 
age do not have time for a single full oscillation before 
optical pumping transfers them to a neighboring lattice 
site. 



TABLE II. Laser parameters used in the simulations and 
the corresponding lattice properties: Detuning 8, Rabi fre- 
quency Q, lattice modulation depth Uo, semiclassical average 
number of oscillations in a lattice site N osc — Q osc t p , and 
saturation parameter so. We use the semiclassical average 
oscillation frequency Q osc discussed in fllq]. Units are given 
in parenthesis and the simulations are labeled by the element 
and the lattice depth. 



O(r) 5{T) U (E r ) N osc s 



Simulation label 



1.2 


-3.0 


374 


0.93 


0.08 


Cs374 


1.5 


-3.0 


584 


0.74 


0.12 


Cs584 


2.5 


-3.0 


1621 


0.47 


0.33 


Csl621 


1.8 


-3.0 


560 


0.76 


0.18 


Rb560 


2.8 


-3.0 


339 


0.98 


0.42 


Na339 


3.5 


-3.0 


530 


0.78 


0.66 


Na530 



Steady state properties of the atomic cloud in the lat- 
tice are characterized e.g. by the average kinetic energy 
per atom, the spatial probability distribution or the mo- 
mentum probability distribution. These are results which 
we obtain from our simulations. We keep the detuning 
fixed (S = — 3r) and vary the Rabi frequency fi, which 
gives various values for the optical potential modulation 
depth 



U = --Ms , 
where Sq is the saturation parameter given by 

n 2 /2 
s ° ~ <5 2 + r 2 /4' 

The spatially modulated optical potentials are 

[/_ = Uq sin 2 (k r z), 
Uq cos 2 (krz), 



U 4 



(2) 



(3) 



(4) 



for ground states m g = —1/2 and m g = +1/2 respec- 
tively | ft8| . The parameters used in our simulations along 
with relevant lattice properties are summarized in Ta- 
ble |n} 

Collisions and radiative heating increase the relative 
velocity between the atoms [|]. This heats up the atoms 
and it is possible for a colliding pair to escape from the 
lattice. One can calculate semiclassically the critical mo- 
mentum p s c c giving the point in momentum (p) space 
where t he co oling force has its maximum value H. In 
Section VI F below, we discuss for which values of the 



momentum p c we may neglect energetic histories and con- 
sider the corresponding atoms lost from the lattice. 



III. BINARY INTERACTIONS BETWEEN COLD 
ATOMS 

In this Section we give a simple description of collisions 
between two cold atoms in the presence of near resonant 
light M and discuss the background of this phenomenon 
to occur in an optical lattice. 
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FIG. 4. Radiative heating of colliding atoms. The quasi- 
molecule is excited at the Condon point r c and accelerated 
on the upper level before spontaneous decay terminates the 
process. 



A. Binary interactions 

In this study, we consider atomic gases with an occu- 
pation density of p = 25 %, i.e., every fourth lattice 
site is occupied and the average distance between two 
atoms is z a — A. For Cs this corresponds to a density of 
1.62 x 10 12 atoms/cm 3 . This atomic gas density is low 
enough that collisions can be treated as binary processes: 
i) The collision range is an order of magnitude smaller 
than z a , and ii) for a Cs atomic mass, atoms with typical 
maximum velocities produced in our simulations need to 
have evolved over a time larger than 75 T -1 to travel a 
distance z a and they scatter a large enough number of 
photons that there is negligible memory effects between 
two collision events. Thus, the binary collision picture is 
justified in our calculations. 

Let us consider two atoms with a temperature around 
or below the Doppler temperature Tjj. If such two atoms 
collide they form a quasimolecule which a near resonant 
light may excite when the atoms approach each other. 
This occurs at an internuclear distance called the Con- 
don point (r c ) where the excited state electronic molec- 
ular potential becomes resonant with the ground state 
potential as displayed in Fig. |j. 

We neglect the hyperfine structure of the atoms and do 
not consider the ground state hyperfine structure chang- 
ing collisions, but concentrate on the effects based on 
radiative heating and escape of the colliding pair. In 
this process, the resonant excitation of a quasimolecule 
terminates in spontaneous decay and the colliding pair 
of atoms gains kinetic energy due to acceleration on an 
attractive molecular excited state before decay occurs, 
see Fig. [|. In principle it is also possible to loose atoms 
from the trap via fine structure changing collisions [0. 
In this work, the loss fraction of fine structure changing 



collisions is assumed to be negligible compared to the 
radiative escape mechanism. The small detuning of the 
lasers makes r c very large, and the probability for sur- 
viving to the small internuclear distances required for the 
fine structure change is rather small. Furthermore, the 
energy increase caused by this mechanism is large and 
leads mainly to loss of atoms, contributing to heating 
only via secondary collisions. These secondary collisions 
are still rather rare in the low-density gas samples of laser 
cooled atoms. 

An essential ingredient in the kinetic energy increas- 
ing collision process is the excitation of the large fraction 
of the population into the attractive excited state. The 
excited state population fraction in turn depends on the 
relative velocity between the interacting atoms when they 
reach r c for the attractive molecular states. For the colli- 
sions in the lattice, the relative velocity in turn depends 
on the optical lattice modulation depth. The deeper the 
lattice, the higher the relative velocity of the atoms when 
they end up in the same lattice site and collide. We con- 
sider here lattices with modulation depths in the range 
339E r < Uo < 1621EV, where E r is the recoil energy, 
cf. Table III. Thus, the relative velocities before a colli- 



sion remain low, which keeps the excitation probability 
large. The small detuning keeps the excitation proba- 
bility large also since the excited state slope decreases 
with the detuning; this increases the interaction time for 
moving atoms at the vicinity of the resonance point r c . 



B. Collisions in lattices 

Theoretical and experimental collision studies in 
MOTs show that the atomic cloud is heated by the radia- 
tive mechanism described above. Atoms may also escape 
from a MOT by this mechanism [§,0. Thus these colli- 
sions set a limit for atomic densities and temperatures of 
the cloud in a MOT when the density is increased so that 
binary interactions begin to have a clear effect. Typical 
densities achieved in MOTs are around 10 11 atoms/cm 3 
and temperatures around or below the Doppler cooling 
limit [g. 

Similar effects are expected in an optical lattice when 
the occupation density of the lattice increases. What is 
not directly expected is that there is a parameter region 
where a possible cooling process in a dense lattice occurs 
due to collisions. This is due to the fact that the collid- 
ing pair of atoms carry more kinetic energy than localized 
ones and during a collision they almost always gain suffi- 
cient energy to escape from the lattice. Thus those atoms 
that have not collided inelastically and remain in a lat- 
tice, have less kinetic energy per atom. Moreover, they 
can also thermalize via elastic ground-ground collisions. 
In our study, we neglect the rescattering of photons and 
consequently the total pattern of cooling and heating is 
not studied. Here we only consider the effects of colli- 
sions. The complete problem is simply not computation- 



ally tractable within our framework. 

Atomic interactions in lattices are usually modelled as- 
suming fixed positions for both atoms and calculating 
how the atomic energy levels are shifted by the interac- 
tion JlQJ [l3[ . These static models ignore the dynamical 
nature of the collision processes described here. When 
allowing the atoms to move the problem becomes compli- 
cated and computationally extremely tedious. To make 
numerically feasible calculations, we have fixed one atom 
and allow the other one to move freely, as described fur- 



ther in Sec. VI B 



IV. 



ATOMIC BASIS FORMULATION AND 
DIPOLE-DIPOLE INTERACTION 



In this section we describe the two-atom product state 
basis [Q and the dipole-dipolc interaction (DDI) be- 
tween two atoms in our one-dimensional (ID) study. 



A. Atomic basis formulation 

We do not use the adiabatic elimination of the excited 
states, which is typically employed in order to simplify 
the equations for atomic motion po| . By keeping the 
excited states in the calculation we are able to account 
for the dynamical nature of atomic interactions and the 
radiative heating/escape mechanism. 

In general the product state basis vectors are: 



\jimi)i\J2m 2 )2, 



(5) 



where j\ and ji denote the ground or excited state (in 
our case g for the ground state 2 Si/2, e for the 2 P3/2 
excited state) and mi, m2 denote the quantum number 
for the component of j along the quantization axis z for 
atom 1 and 2 respectively. The total number of states is 
6 x 6 = 36. 

We have t o fix the position of one atom, as described 
in Sec. VI B. If the position of atom 1 is fixed, the binary 



system wave function depends now only on the position 
of the moving atom 2 



\Tp{Z2,t)) 



E 



V , i 1 ,'m, 1 ( z 2,i)|jimi)i|j 2 m 2 )2- (6) 



31 J2,m 1 ,m 2 



The atomic spatial dimensionality of the problem is re- 
duced from two to one. The relative coordinate z between 
atoms is now z — z 2 — zt where Zf is the position of the 



VI B 



fixed atom, see Sec. 

In the atomic product state basis 
Hamiltonian is 



H S = H 1 +H 2 + V d 



I {)• 



]19||, our system 



(7) 



Here, Vdip includes the interaction between the atoms and 
Hi = H\®^2 and i/ 2 = li<g>-ff 2 , where the operators l a 



are unity operators in atom a subspace, and the single 
atom Hamiltonian for atom a (a — 1,2) is, after the 
Rotating Wave Approximation (RWA), 



S a = ^-nsp e , a + v a . 



(8) 



Here, P eyCl — i^ m= _ Z j 2 



&m)a a(e m |, and the interaction 



between a single atom a and the field is 

V a = -i— =sin(/cz Q ){|e 3 /2)a 0(01/2! 

+— | |ei/ 2 )a a(9-l/2\} 

+ —^=COs(kz a ) {|e_ 3 / 2 )a a (5-1/2 1 

+ ^ l e -l/2)a a (#1/2 I } + h.C, 

where z a is the position operator of atom a. 



(9) 



B. Resonant dipole— dipole interaction 

In order to get the DDI potential, Vn p , in Eq. (ffl), 
we have calculated the Master Equation for the atom 
and laser field in question, and identified the DDI. Our 
approach follows the lines of Appendix A in Ref. pl| . 
As it is beyond the scope of this paper to go through 
the derivation of the DDI potential in detail, we shall 
refer to equation (Ax) in Ref. Q as Eq. (LMAx). We 
identify Vdi p as the terms similar to An and A22 with 
{n u + 1) = 1 in Eq. (LMA21). 

First, it is convenient to write the non-interacting sys- 
tem Hamiltonian H\ + iJ 2 in a basis of center-of-mass 
and relative coordinates: 



P = Pl+P2, P 



P2-P1 



(10) 



With these coordinates, the interaction potential with 
the laser field, V = V\ <8> I2 + H-i <8> V2, reads 

v - m ■ 

V = —1 — = sin 



V2 
Ml 



(kZ) cos (fc|) (Sl t+ <g> l a + Hi <8> 5+ >+ ) 
+i-i= cos(fcZ) sin (k^\ (S+ j+ ® 1 2 - li ® 5+,+) 
+i—7= cos(fcZ) cos (k^\ (S\_ (g) 1 2 + li ® 5+ _) 

+i^= sin(fcZ) sin (fc-J (5}__ ® 1 2 - li ® 5+ _) 

+/i.c, (11) 

where Z and z are the center-of-mass and relative coor- 
dinates along the z-axis and 



m=l/2 

+ ,9 = Z-^i CG m \e m +q)a a{gm\- 
m=-l/2 



(12) 



Here CG^ are the appropriate Clebsch-Gordan coeffi- 
cients and q is the polarization label in the spherical ba- 
sis. We rewrite similarly the interaction with the vac- 
uum electromagnetic field in terms of the relative coor- 
dinates ( J10| ) . The DDI terms are identified after we have 
considered the damping part of p [cf. Eq. (LMA17)] in 
the derivation of the Master Equation for our two-atom 
system. 

Following pl| , we note that the DDI potential is found 
as 



3 1 Z" 00 / ui 

Vdi P = -„ftr- / duj[ — 

8 7T J \UJ 
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^UJ — UJq 



-S+qS-q 



P 2 (cos9 r )[--(S ++ S 



-5_i_n5. 



+0'- , -0 



1 



-S+-S-) 

,, „, P2(cos8 r )cos(t) r x (13) 

(— S++S-0 + S+oS — — S+oS — h + 5^ — S-o) 
+ -P 2 2 (cos(9 r )cos20 r . {S++S—+S+S-+) }, 

where P(x) is Cauchy's principal value, ji are spherical 
Bessel functions of the first kind, P2 is Legendre poly- 
nomial, and P^ are associated Legendre functions. The 
angles r and <f> r are the angles of the relative coordinate 
r in the spherical basis. We have also introduced the 
operators 



S+qS-q' — (S +q S_ q , 



S+,qS-,q') 



(14) 



Thus, we need to calculate integrals of the type 



where S" „ 



1 f°° ( lo \ 3 
Mvor) = - / du[ — )P 
k Jo K^oJ 



uj — UJq 



Ji: 



H) ' ^ 



with qo = Wo/c. We may change the lower limit in the 
integral to —00, enabling us to calculate the integral by 
contour integration [B2j . The results are 



Zofaor) 



cos <7o r 



Zifaor) = - 



Qor 
cos qor 



(16) 



sm q r cos q$r 



Qor \ (qor) 2 (q r) 3 

and the three-dimensional DDI potential is 

V dip = -hr (i^^(l - 2P 2 (cos0 r ))x 
o {J qor 
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FIG. 5. The shifted ground state and the attractive ex- 
cited state [labeled by Hund's case (c) notation] molecular 
potentials of Cs for 5 — — 3.0I\ 
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If the two atoms are positioned on the z-axis, the DDI 
potential reduces to the one-dimensional potential 
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By diagonalizing Vdi P it is possible to obtain the molec- 
ular potentials shown in Fig. a. One also notes that the 
DDI induces the 7r-polarization couplings which the laser 
fields do not do here. 



V. MONTE CARLO WAVE FUNCTION METHOD 

In this section we describe briefly the main features 
of the Monte Carlo wave function (MCWF) method Q 
which was developed for problems in quantum optics and 
discuss the implementation of the method to solve the 
cold collision problem in optical lattices. 



A. Basic Monte Carlo method 



Various types of Monte Carlo (MC) methods [I4|-|l6|] 
have been developed for problems where a direct ana- 
lytical or numerical quantum mechanical solution of the 
density matrix Master Equation is very difficult or im- 
possible due to the complexity of the problem. Complex- 
ity usually arises because of the coupling of the system 



studied to a reservoir with a large number of degrees of 
freedom and also because of a large number of elements 
in the system density matrix. Problems of this kind are 
common in laser cooling of neutral atoms. Various types 
of quantum approaches are possible in 2D systems |23j ] 
but in 3D a full quantum treatment of laser cooling of 
atoms has only been given in terms of the Monte Carlo 
method ||. 

The core idea of MCWF method is the generation of a 
large number of single wave function histories including 
stochastic quantum jumps of the system studied. So- 
lutions for the steady state density matrix and system 
properties can then be calculated as ensemble averages 
of single histories. 

To generate single histories of the system wave function 
\i[>), one solves the time-dependent Schrodinger equation 



Here the non-Hcrmitian Hamiltonian H is 
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where Hs is the system Hamiltonian, Eq. (R) in our case, 
and the non-Hermitian part Hdec includes the decay 
part. Hdec is constructed from appropriate jump oper- 
ators Cj corresponding to a decay channel j and to the 
detection scheme of the system. The general form of the 
non-Hcrmitian part reads 



H 



DEC 
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(21) 



During a time evolution step dt the norm of the wave 
function may shrink due to Hdec and the amount of 
shrinking gives the probability of a quantum jump to oc- 
cur during the short interval St. Based on a random num- 
ber one then decides whether a quantum jump occurred 
or not. Before the next time step is taken, the wave 
function of the system is renormalized. In the case that 
a jump occurs, one performs a rearrangement of the wave 
function components according to the jump operator Cj, 
corresponding to decay channel j, before renormalization 

of |V>- 



B. Time evolution 

A natural termination for a simulation occurs if a 
steady state appears. One must note that the time to 
reach steady state varies even for the same system stud- 
ied when the laser parameters are changed. Therefore 
one has to be careful to have long evolution times to 
make sure that the steady state is reached and ensemble 
averaging can be done in a reliable way. 

We solve the time-dependent Schrodinger equa- 
tion Eq. Jl9| ) by the split operator-Fourier transform 
method |25||. Formally solving Eq. (|l^) over St gives 



m + st)) = um )) 

where the time evolution operator U reads 



U = exp 



iHSt\ 



(22) 
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We split the time evolution operator U including the 
Hamiltonian H of Eq. ( f20| ) into three parts as H = Hy + 
Hk + H]j. When H is in matrix form, Hy has an off- 
diagonal part accounting for the atom-field coupling and 
the interaction between atoms, Hk is the diagonal kinetic 
part and Ho includes the non-kinetic diagonal part, i.e., 
decay and detuning. 

For non-commuting operators A and B we can write 
to second order accuracy 25 



exp (A + B) ~ exp {A/2) exp (B) exp(A/2). (24) 

As we take many consecutive time steps during the evo- 
lution, we finally approximate the wave function at time 

to + nSt by 



\ip(z,to+ndt)) 



n 



U V U D /2 U K U D /2 



hK-Mo))- (25) 



Here, UD = e*p(—'iH D 5t/h) and XJk = 3 : ~ 1 exp(— iSt^jj) 
where T and T~ x denote the Fourier and inverse Fourier 
transforms. Finally, Uy can be written as Uy — 
S e~scp(D)S~ 1 where S contains eigenvectors and D eigen- 
values of Hy. Uy corresponds now to a change of basis, 
multiplication by exponentials of eigenvalues and change 
of basis back to the product state basis. The above form 
for the temporal evolution of \?p) is straightforward to 
implement and fast on a computer, e.g., 20% faster than 
the Crank-Nicholson method. 



C. Decay channels 

A single atom has six different ways to spontaneously 
emit a photon so the total number of decay channels is 
12 for two atoms (channels f — 6 for atom 1, 7—12 for 
atom 2). For each decay channel, j, the jump probability 
is given by 



Pi = 5mc\c 3 



(26) 



where the jump operators Cj are constructed from single 
particle jump operators. In the single particle subspace 
for atom a and decay channel j we have 



C Q = CGjVT \g a m a ) a a {e a m a \, 



(27) 



where e a m a labels the excited level from which, and 
g a m a the ground level to which jump occurs. Exten- 



sion to the product state basis is simple M 
1, Cj = Cj (g) 1 2 , and for atom 2, Cj = l x ® 



For atom 



For example, if we denote the jump of atom f from 
l e -i/2)i to |.g-i/2)i as channel 2, the jump operator in 
the product state basis for this jump is 

C 2 = y / 2/3Vf{\g^ 1/2 )i 13-1/2)2 i(e-i/ 2 | 2(9-1/2! 
+1.9-1/2)1 1.9+1/2)2 i(e_i/ 2 | 2(9+1/2! 
+l.9-i/2)i |e_ 3 /2)2 i(e_i/ 2 | 2(e_ 3 / 2 | 
+ l.9-i/2)i |e-i/2>2 i(e_i/ 2 | 2(e_i/ 2 | 
+|9-i/2)i |e+i/ 2 ) 2 i(e_i/ 2 | 2<e + i/ 2 | 
+ |9-i/2)i |e+ 3 / 2 )2 i(e_i/ 2 | 2(e +3 / 2 |} (28) 

and the corresponding jump probability for channel 2 is 
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We neglect here the case where both atoms jump and 
two photons are detected simultaneously. The proba- 
bility for a single atom jump during St is <C 1 so the 
joint jump probability is negligible compared to the single 
atom jump probability. In principle it would be possible 
in simulations to take into account joint jumps but this 
unnecessarily complicates the jump procedure. After ap- 
plying the jump operator Cj, the wave function is still in 
a superposition state, but it has collapsed onto product 
state basis vectors, leaving only one ground state level 
component of the jumped atom populated. 



D. Ensemble averaging 

We calculate the results as an ensemble average of sin- 
gle history time averages in the steady state time do- 
main [ fL6[ . This averaging method requires a smaller 
number of histories calculated to achieve reasonable er- 
ror bars than a simple ensemble averaging at a single 
steady state time point. In our study, extra complica- 
tions arise because the number of collision processes in 
the whole ensemble also comes into play. Atoms do not 
end up in the same lattice site, producing collisions, in 
all the calculated histories. Atomic hopping between lat- 
tice sites is a stochastic process and only in a fraction of 
the total number of histories, collision processes occur. 
We need a sufficiently large ensemble to produce enough 
collision events to have reliable results. This is why we 
have a much larger ensemble size, 96— 128 members, than 
e.g. used in 3D laser cooling Monte Carlo simulations us- 
ing the same ensemble averaging method [|24j . 

In order to be able to use this averaging method, we 
need to be sure that time averages of single histories have 
been calculated in the steady state time domain. The 
simulation times used are displayed in Table IV . We have 



<?• 



carefully checked from the time evolution of the kinetic 
energy that the simulation time was long enough to reach 
well into the steady state time domain. 



TABLE III. Characteristic units. Distances are given in 
nm, momenta in 10 -28 kgm/s, time in ns, and energy in 1CP 30 
J. 



Quantity 


Characteristic unit 


Cs 


Rb 


Na 


distance 


A = 1/kr 


136 


124 


94 


momentum 


T)f /ZfCy 


7.77 


8.49 


11.25 


time 


r- 1 


193 


167 


99 


energy 


E r = {h 2 kl)/2m 


1.37 


2.55 


16.57 



TABLE IV. Simulation times: Total time for collision sim- 
ulations T, ensemble averaging time T aV e, time step size St 
and maximum momentum |p| moa; given by 8t for the numer- 
ics to remain reliable. 



Simulation 


T(r _1 ) 


lave\i- Sq J 


^(r- 1 ) 


\P\ max \Pr) 


Cs374 


1600 


78-125 


0.2 


110 


Cs584 


1600 


97-194 


0.2 


110 


Csl621 


760 


178-256 


0.1 


155 


Rb560 


1600 


140-280 


0.2 


89 


Na339 


470 


99-198 


0.05 


90 


Na530 


470 


155-311 


0.05 


90 



VI. SIMULATION SCHEME 

In this section we present the charasteristic units used 
in our calculations and discuss the various criteria which 
set the numerical limits for the simulations. The approx- 
imations used and the numerical details related to the 
wave packet initial conditions and dynamics are also pre- 
sented. 



A. Scaling and discretization of space 

It is useful from a practical point of view to choose 
suitable units and scale the time-dependent Schrodinger 
equation ( |l9|) accordingly. Convenient physical units and 
their numerical values for the th ree appropriate alkali 
metal elements are listed in Table III. In the discussion 



below, we list all quantities and scale equations in units of 
the characteristic quantities displayed in Table III unless 
explicitly stated otherwise. 

As the phase factor exp(— iE5t/e) has to be well de- 
fined, cf. Eq. (|23|) , we obtain a criterion for the maximum 
size of the time step St dictated by the maximum kinetic 
energy since we should fulfill the relation St <C e/p 2 - Here 
e is the energy of the linewidth of the transition in re- 
coil units, as previously displayed in Table QL Collisions 
increase the atomic kinetic energies which makes the cri- 
terion for St numerically more strict for the two-atom 
case, compared to the one-atom Sisyphus cooling simu- 
lation (cf. Fig. ||). We give in Table |y| the values of St for 
various simulations and the maximum momentum |p| maa; 
for the numerics to remain reliable. The total simulation 
times are 125 — 311 in units of 1/(Tsq). These depend on 



the properties of the alkali metal elements and the laser 
parameters. 

For the numerical simulations, one has to discretize the 
position and momentum spaces, and the resolution has 
to be fine enough to ensure valid results. We have used 
8192 grid points when the length of the entire spatial grid 
is L z = 5 A ~ 31. 4A. This gives the step sizes in position 
and momentum spaces of Sz ~ 0.0038 and Sp = 0.2. 
The width of the Gaussian wave packet at the begin- 
ning of the simulation is Azq = 0.02A ~ 0.1257A giving 
Azq/Sz ~ 33 and in momentum space Ap /Sp ~ 20, 
ensuring sufficiently fine resolution. 

The inverse space (here momentum space) has reflect- 
ing boundary conditions when using the FFT method. 
Thus the size of the momentum space has to be large 
enough to avoid the reflection of high kinetic energy 
atoms at the edges of the momentum space. This re- 
quires special attention when considering the interaction 
simulations where the kinetic energies of the atoms in- 
crease due to the inelastic collisions. 

The momentum space grid has a total size of L^ = 
2ir/Sz = 1638 so that the atomic momenta may have 
values \p\ < 819. The depths of the lattices in our sim- 
ulations are such that atoms localized at a lattice site 
have momenta \p\ < 50. The momenta increase when 
the atoms wander around in the lattice, especially due to 
the inelastic collisions. The probability of gaining a suf- 
ficient momentum to reach the edge of our momentum 
space grid in a single collision event is now negligible. 
On the other hand many consecutive collisions do not 
shift the population for large p and the reflection effect 
is avoided. This is due to the fact that the increasing 
relative velocity between the atoms reduce the excitation 
probability and increases in momentum terminate before 
the edges of the p-space are reached. 



B. Position fixing of one atom 

We simulate the behaviour of a 36 level dissipative 
quantunr system with a position-dependent coupling to 
the laser field and a position-dependent coupling between 
two atoms. This requires large computational resources. 
With the current computer capacity, it is not possible to 
simulate the situation where both atoms are allowed to 
move freely. Instead we have to fix one atom spatially and 
let only the other atom move. This reduces the dimen- 
sionality of the problem to one since the relative position 
of the atoms with respect to the laser field is now fixed. 
This also means that an inelastic interaction process will 
not change the kinetic energy for both atoms, but we use 
the relative kinetic energy as an estimate for the kinetic 
energy change per atom. 

In our previous study |9| the position of the fixed atom 
was kept constant but here we relax this condition. The 
position Zf of the fixed atom is now selected randomly 
in the interval \zf\ < 0.1 25 A for each ensemble member. 



This range covers all the interesting physics as an atom 
fixed outside the current range would be rapidly opti- 
cally pumped to the opposite lattice well and the situa- 
tion would correspond to that with the above-mentioned 
range of Zf . 

The change of Zf also moves the Condon point r c with 
respect to the lattice and now r c may be located any- 
where between the lattice well and peak. This is an im- 
portant point for making our model more realistic: Since 
the kinetic energy of the atom changes when it moves 
in the periodic optical potential, the relative velocity be- 
tween the atoms and thus the excitation probability to 
the attractive molecular state at r c depends on its posi- 
tion in the lattice. When r c is located at the peak of the 
optical potential, the atom has to move up the potential 
hill to reach it. The relative velocity between the two 
atoms is now less and the excitation probability higher 
than in the case where the location of r c is at the bottom 
of the optical potential. 



steady state after cooling in the lattice does not neces- 
sarily correspond to a Maxwell-Boltzmann distribution 
of velocities but a clear steady state corresponding to the 
lattice properties is still reached jl8[ . 

It should be noted that although we include by default 
the recoil effects by absorption and stimulated emission, 
we cannot in our quantum approach take proper account 
of the Doppler shift which is the basis of the semiclassical 
description of Doppler cooling |^7j . So in practice we ne- 
glect the Doppler cooling mechanism. Adding the recoil 
from spontaneous emission will not change this fact, nor 
the simulation results. In any case, the role of Doppler 
cooling is negligible when compared to the Sisyphus cool- 
ing force for the velocities of atoms localized in the lat- 
tice j|. It can, however, have a role in the recapture of 
the hotter atoms heated by collisions, so in that sense 
our model is limited. 



D. Occupation density of lattice 



C. Initial wave packet 

At the beginning of the simulation, the wave packet is 
placed in a randomly chosen ground state sublevel. The 
initial Gaussian packet has a full spatial width of 0.02A. 
Thus the initial position of the spatially relatively narrow 
wave packet in the lattice is random but completely out 
of the range of the molecular resonance. 

Each wave packet has a randomly selected mean ini- 
tial velocity given by the Maxwell-Boltzmann distribu- 
tion corresponding to the selected initial temperature of 
the atomic cloud. We emphasize that the momentum 
space width of the initial wave packet has no associa- 
tion with the thermal distribution, as it is merely needed 
to satisfy the Heisenberg uncertainty relation for a spa- 
tially localized initial state. As stated above, the con- 
nection between the wave packet and the temperature 
takes place via the mean momentum of the wave packet. 
By selecting this mean momentum randomly for each en- 
semble member but weighting the occurrences with the 
Maxwell-Boltzmann distribution, we create within the 
Monte Carlo ensemble another ensemble of possible ini- 
tial collision velocities. This is the wave packet version 
of the standard collisional energy average 1 26 1 . 

As mentioned above only the ensemble averaged mo- 
mentum probability distribution has a relation to tem- 
perature. This initial distribution gets narrower when 
the system evolves and the simulation progresses corre- 
sponding to cooling of the whole atomic cloud. Moreover, 
it must be stressed that the steady state reached does not 
depend on the initial widths of the single wave packets 
nor on the initial temperature as long as the atoms are in 
the reach of Sisyphus cooling. The simulation times get 
longer when the initial temperature is increased but we 
want to take into account the effect of collisions on the 
cooling dynamics in a realistic way. We also note that the 



We have performed all the simulations presented in 
this paper for an occupation density of p = 25% in one 
dimension, whereas in our earlier work |9[ , we presented 
also results for other one-dimensional densities in Cs lat- 
tices. These previous simulations showed that an occu- 
pation density of 25% is sufficiently high for interesting 
effects to appear, namely an evaporative cooling process 
which works for at least some parameters of the laser 
field. The occupation density used in this paper is also 
nearly the largest density we can use when the simula- 
tions are done in the way presented here. The purpose 
of the paper is to further explore the parameter space 
and to extend our simulations to other atomic masses in 
addition to describing the details of our simulation ap- 
proach. 

For p — 25% the available spatial length for the mov- 
ing atom should be equal to A, corresponding to the aver- 
age distance between the atoms. But decreasing the spa- 
tial size increases the step size in the momentum space 
since 5p = 2tt/L z . So to have a sufficiently fine resolution 
in momentum space and still keep p — 25%, we choose 
L z = 5A and set an elastic repulsive potential barrier 
such that the allowed spatial length is A. The forbidden 
spatial region thus makes the numerics work properly 
without altering the physics. Of course the forbidden 
region does not affect the momentum space. 

Every fourth lattice is occupied when p = 25% or in 
other words there is 1 atom per wavelength A. This cor- 
responds to a situation where the fixed atom sits e.g. at 
z = and the repulsive elastic potential barrier for the 
moving atom is set at z = A. Then two consecutive col- 
lisions are described when the moving atom travels from 
the first collision region to the repulsive barrier, turns 
back, and collides again. Memory effects from previous 
collisions are rapidly removed due to decoherence, as also 
discussed in Sec. III. This is why we can say that the 
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present model describes collisions in general in a lattice 
and not only between the same two atoms. 

The dimensionality of the problem and the position 
fixing of one atom causes subtleties related to the occu- 
pation density: The moving atom travels on average a 
distance of 1A for the first collision. After this event it 
has to travel a distance 2A (from Zf to Zf + 1 and back) 
to collide again. The probability is high for a large ki- 
netic energy increase during the first collision. Thus the 
first collision has the dominant effect on the kinetic en- 
ergy scale relevant to the lattice dynamics. This is why 
we rather use here p = 25% which corresponds to the 
average collision distance of 1A. 

It is not trivial to connect the one-dimensional occu- 
pation density to the two- or three-dimensional density, 
as that will depend on the particular form of the lat- 
tice. It suffices to note, though, that a certain occupa- 
tion density in one-dimension will typically correspond to 
a lower value in higher dimensions, so that e.g. in three 
dimensions we expect to see the effects studied here for 
three-dimensional occupation densities less than 25%. 



E. Interaction at short range 

The DDI , Eq. (|l8|), becomes singular at short range. 
The singularity in Hy is simply removed by replacing r 
with r+r ff and choosing r ff = 1CT 8 . When constructing 
the time evolution operator U, we diagonalize Hy and 
the DDI part in Hy produces the eigenvalue manifolds 
corresponding to the attractive and repulsive molecular 
potentials. 

We replace the position dependencies of the attractive 

Eq. ©, by 
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states which are the same as in V£? ls 
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(r b + r off ) n/b 



, (n= 1,2,3) 



(30) 



in a manner similar to what was done in Ref . |7j . Table |V] 
gives the used values of 6, and we show the potentials in 
Fig. |. 

The main reason for this "flattening" of the attrac- 
tive state potentials is that considering the numerics, we 
have an upper limit to momentum. Thus we need to set a 
maximum momentum which can be reached in our simu- 
lations by acceleration, but which can still be treated reli- 
ably numerically in our integration grid, and is neverthe- 
less large enough to correspond to a clear loss process. By 
selecting different values for b for each molecular poten- 
tial we take into account the individual characteristics of 
the different attractive states and of the atomic elements. 
It should be noted that by the time the atoms reach the 
artificially modified part of the attractive potentials, they 
move fast enough to make decay unlikely before they are 
reflected and move again to the region where the modifi- 
cation does not affect the potentials Thus the flattening 
of the potentials does not increase the time the atoms 
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z (V2ji) 
FIG. 6. The short range attractive excited molecular po- 
tentials for Cs. Repulsion of exponential form has been added 
and the deep part of the potentials flattened to allow reliable 
numerical treatment of momentum (see text). 

TABLE V. Values of b. Numerical values which are used 
for the various attractive molecular states and atomic species. 
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spend inside the modified part of the potential, i.e., this 
does not enhance radiative decay artificially. 

We concentrate on the radiative heating which occurs 
because of the strong decay in the vicinity of r c . In this 
region the treatment of the singularity in the DDI and 
in the attractive molecular states does not yet affect the 
potentials. For example, the removal of the singularity 
causes a change in the value of the Cs 0+ potential of 
0.7% at the position r = 0.50 when r c ~ 0.83 for the 
detuning used, S = —3T. 

The atoms repel each other at the very short range 
when their electron clouds begin to overlap. We do not 
have to consider the details of the short range repulsion. 
Thus the short range repulsion is simply produced by 
adding the exponential term a exp(— f3r) e to the eigen- 
values of Hy . For flat states (states other than attractive 
or repulsive excited statesVa = 30, (3 = 20 for Cs and Rb, 
a = 100, /3 = 20 for Na p8[ . For the attractive excited 
state eigenvalue manifolds a — 25, j3 = 15 for Cs and 
Rb, a — 90, j3 = 15 for Na. Values of a and (3 are chosen 
such that they produce high enough repulsion within a 
sufficiently short range but without producing numerical 
difficulties because of the contradictory requirements of 
height and range. 

Finally we emphasize that apart from flattening the 
potentials, we perform our calculations in the atom- 
atom basis. Thus the molecular states are not used di- 
rectly. Unfortunately our dipole-dipole potential takes 
care only of the force part of the atomic interaction. The 
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r-dependence of the lifetimes of the molecular states are 
ignored, i.e., each molecular state ends up having the con- 
stant atomic linewidth, instead of the retarded linewidth. 
This linewidth arises from the fact that the two atoms 
couple to the same vacuum modes at different locations, 
leading to an e~ lk ' r phase difference term. However, the 
atomic lifetimes differ at maximum only by a factor of 
2 from the atomic one. The main exception is the 2 U 
state, which becomes strongly dipole-forbidden and can 
support strong survival and thus e.g. favor the fine struc- 
ture change loss mechanism over the radiative process. 
But we have typically r c < 0.8A, which means that the 
2 U state is already hard to excite as well (see Ref. |2£]] 
for a detailed discussion). In order to make the quantum 
jump process tractable we use the atom-atom basis, and, 
unfortunately, we can not just transform into the molec- 
ular basis, change the linewidths into retarded ones, and 
then transfer into the atom-atom basis (as we do when 
flattening the potentials). This is because for decay, the 
lifetimes appear in the jump operators in addition to the 
Hamiltonian. 



F. Atoms escaping the lattice 

One needs to define the critical momentum p c to be 
able to calculate the MC ensemble averages and the prop- 
erties of the atoms remaining in the lattice. If \p\ < p c 
atoms are considered to remain/relocalize in the lattice, 
whereas when \p\ > p c they are considered lost from the 
lattice due to a collision or a series of collisions. Semiclas- 
sically, we can calculate the critical p s c c where the cooling 
force has its maximum value for the parameters used [Q . 
The cooling force is, of course, still effective for momenta 
above p s c c . 

Due to the stochastic nature of the jumps, it is not 
possible to say if a given high momentum atom will re- 
localize or if it is lost from the lattice. Assume that an 
atom has a momentum |p| > p s c c . If the following few 
quantum jumps reduce the kinetic energy of the atom, 
depending on the atomic position in the lattice, it has a 
good chance to relocalize in the lattice. In the opposite 
case, where the next few jumps increase the kinetic en- 
ergy of the atom, corresponding to jumps from the vicin- 
ity of the bottom of the potential well to the vicinity of 
the top of the well, the atom has less probability to relo- 
calize into the lattice. This means that for two different 
MC histories with the same initial value of \p\ > p^ c one 
atom may escape from the lattice whereas the other one 
may relocalize. Thus it is not possible to define p c in a 
way that all atoms below p c always relocalize while when 
|p| > Pc they escape. 

When we calculate the kinetic energy per atom stay- 
ing in the lattice, we need a criterion for neglecting those 
MC histories in the ensemble averaging that correspond 
to atoms lost from the lattice. To solve the problem, we 
have calculated the kinetic energy per atom by using var- 
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t (r- 1 ) 
FIG. 7. Kinetic energy time evolution indicating a too 
large choice for p c . The steady state is not reached since col- 
lisions increase the kinetic energy and the collided atoms are 
out of the recapture range, and thus escape from the lattice. 
(Cs584, Pc = 60). 
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t(r- 1 ) 
FIG. 8. Kinetic energy time evolution indicating a steady 
state. Atoms above the recapture range are now omitted with 
the proper choice for p c . (Cs584, p c = 40). 



ious values of p c . Since there is an increase in the average 
kinetic energy as a function of time when the value of p c 
used is too large, we may check from the time evolution 
of the kinetic energy that our choice for p c is the proper 
one when we want to calculate the average kinetic energy 
per atom in the lattice. This is because more collisions 
occur as the system evolves in time and if the gain in 
kinetic energy is too large for the atoms to relocalize in 
the lattice, the kinetic energy increases as a function of 
time and no steady state is reached, as demonstrated in 
Fig. ffl. Whereas when we use an appropriate value for 
p c , atoms still relocalize in the lattice and the kinetic en- 
ergy exhibits a steady state behaviour, cf. Fig. pi It is at 
the transition point between these two different types of 
behaviours of the kinetic energy that we should choose 
the correct value for p c . 

Consequently, the atoms of a kinetic energy exceeding 
the limit given by p c are neglected when we perform the 
ensemble averaging to find the result for the kinetic en- 
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ergy per atom remaining in the lattice. It is important 
to note that the main result related to the narrowing of 
the momentum probability distribution due to collisions 
still includes all the calculated histories and is totally 
independent of p c . 

If there is no constant injection of atoms into the lat- 
tice, collisions slowly deplete it. Finally the density is 
sufficiently low that the interactions between atoms are 
negligible and the atomic cloud regains the properties de- 
termined by the laser parameters only. It should thus be 
realized that what we describe here is a temporary cool- 
ing process which is not effective when the density has 
decreased. What we emphasize is the unexpected be- 
haviour of the system in the intermediate regime where 
the effect of collisions is not heating but cooling. This 
does not represent the nature of the complete dynamics 
of the atomic cloud, of course, as there are other mecha- 
nisms, such as the radiation pressure from scattered pho- 
tons, for both heating and cooling, which are not included 
here. 



G. Computational resources 

The numerical simulations are demanding since we are 
dealing with a 36 level quantum system including various 
position dependent couplings and dissipative coupling to 
the environment. We use 32 processors of an SGI Origin 
2000 machine which has 128 MIPS R12000 processors of 1 
GB memory per processor [B9J . The total memory taken 
by a single simulation (fixed S, fi, p and atomic species) 
is 14 GB and generating a single history requires 6 hours 
of CPU time. A simulation of 128 ensemble members 
then requires a total CPU time which is roughly equal to 
one month. The normal clock time is, of course, much 
shorter (roughly 22 hours) since we take advantage of 
powerful parallel processing. 



VII. THE SEMICLASSICAL APPROACH 

In this section we describe the semiclassical approach 
to calculate the excitation and survival probabilities on 
the molecular excited states of our two-atom system. 



A. Landau— Zener formula and classical path 
approximation 



TABLE VI. Semiclassical probabilities to gain kinetic en- 
ergies on various attractive excited state molecular levels, see 
Fig. ra, for our Cs584 simulation. The probabilities are cal- 
culated for ptot = Per + Ap = 40 and p cr — 24, this value of 
p cr corresponds to the lattice depth. The wave packet spends 
a time t > t e on the excited state, Plz is the Landau-Zener 
excitation probability and P a gives the survival probability. 
Ptot — PaPhz is the total semiclassical probability for the 
process p > p to t = 40 to occur. 



potential 


M^ 1 ) 


Plz 


Pa 


Ptot 


1„ 


0.21 


0.71 


0.81 


0.57 


2u,0g 


0.40 


0.88 


0.67 


0.59 


1 S 


0.48 


0.91 


0.62 


0.57 


0+ 


0.49 


0.92 


0.61 


0.56 



r-dependence of the excited state (C^/r 3 ) Q. The basic 
idea here is that the short resonance region is approxi- 
mated to consist of two spatially linearly behaving states 
and each component of the wave packet arriving at the 
resonance region is independently excited. A more de- 
tailed description can be found in Ref. Q. 

One can then calculate the time it takes to reach a 
point r on the excited state by using the classical path 
approximation 



t = t(r) = - / dr' 



l(& + %- M 



m V 2m r' : 



-1/2 



(32) 



where p cr is the momentum at the Condon point r c . 
There is a direct correspondence between the reached 
point r on the excited state and the energy gain while ac- 
celerating on the attractive excited molecular state. By 
using Eq. (32) one can calculate classically how long it 



takes to reach a point r corresponding to a given increase 
in kinetic energy or momentum due to the acceleration 
on the attractive excited state. 

It is now easy to numerically calculate the kinetic en- 
ergy increase due to the collisions if the wave packet stays 
on the attractive excited state for a time corresponding 
to the natural decay time T -1 . When the exponential 
decay from the excited state is also taken into account, 
the probabilities for various kinetic energy gains due to 
collisions as a function of relative velocity of the colliding 
atoms at r c may also be calculated M]. 



B. Post collision momentum in the lattice 



One can calculate the semiclassical excitation proba- 
bility of the wave packet traveling through the crossing 
region between the two states of the system by using the 
Landau-Zener formula 



Plz = 1 - cxp(-TrA) 



(31) 



where A includes both the coupling between the two 
states and the C3 factor which gives the inverse cubic 



We obtain the values of C3 for the attractive poten- 
tials by fitting near the resonance region the simple ex- 
pression — C3/?* 3 to our molecular potentials obtained by 
diagonalizing the dipole-dipole coupling presented in the 
two-atom basis, Eq. (lq). Figure |] shows an example 
of the total post collision momentum p tot — p cr + Ap 
as a function of p cr when the wave packet spent a time 
t = r _1 on the l u excited state (Cs584). One notices 
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FIG. 9. The total momentum ptot = Per + Ap as a func- 
tion of resonance point momentum p cr for the attractive l u 
state by the semiclassical calculation. The wave packet has 
spent a duration corresponding to F _1 accelerating on the 
excited state before spontaneous decay back to the ground 
state. (Cs584). 
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FIG. 10. Momentum probability distributions for inter- 
acting and non-interacting cases. (Cs374). 
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that Pcr = 24 corresponding to the lattice depth used 
already gives a total momentum of ptot = 67 after a col- 
lision thus pushing the atom to the region in momentum 
space where its probability for rclocalizing back to the 
lattice is small (p s c c — 16.2). This shows in a clear way 
that increases in kinetic energy that are large compared 
to the latttice modulation depth Uq may occur on a time 
scale of r _1 . 

Moreover, when the exponential decay and PlZi 
Eq. (J3lj), are taken into account, one is able to calcu- 
late the probabilities to gain various amounts of kinetic 
energies due to the collision. The total probability P tot 
for the atomic momentum to have at least the value ptot 
after the collision is 



P tnt — r„r\ 



LZ 



(33) 



where P a gives the survival probability on the excited 
state. An example of the results of Ptot are shown in 
Table [Vj. This suggests that the first resonance molecu- 
lar potential 0+ has a dominant role in collisions. For 
Ptot — 40 the probabilities for the various states are 
roughly equal but if the first resonance potential excites 
and accelerates half of the colliding atoms only half of 
them is left for the remaining potentials. 

The simple semi-classical calculation above is not able 
to give quantitative results but it shows that the prob- 
ability to produce an atom of large momentum due to 
a collision is high already when we consider one excited 
level only. This probability increases when we take into 
account that during one collision process, the molecule 
may be excited at four different values of r c related to 
five different attractive states. 
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FIG. 11. Momentum probability distributions for inter- 
acting and non-interacting cases. (Cs584). 



VIII. SIMULATION RESULTS 

The calculated numerical values of kinetic energy per 
atom for various simulations arc shown in Table VII 



and corresponding momentum probability distributions 
in Figs. |H]-[l|. 

Most of the simulations with the selected parameters 
produce a reduced value for the kinetic energy per atom 
when the interactions between the atoms are taken into 



account, see Table VII. Since the inelastic collisions here 



always increase the kinetic energy of the atoms via the 
radiative heating mechanism, our results suggest that the 
consequence of collision almost every time is the escape of 
the colliding high energy atoms from the lattice |3(j] . The 
atoms left in the lattice then have less average energy. 
This is due to the fact that the more energetic atoms 
are favoured to participate in the collision process due to 
their better ability to move between the lattice sites. 

The cooling process is indeed observed when looking 
for the momentum probability distributions includin g al l 
the MC histories for ensemble averaging, see Figs. [lCJ- 
[fq. One can see the slight narrowing of the momentum 
distributions corresponding to the cooling process. The 
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FIG. 12. Momentum probability distributions for inter- 
acting and non-interacting cases. (Csl621). 
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FIG. 14. Momentum probability distributions for inter- 
acting and non-interacting cases. (Na339). 
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FIG. 13. Momentum probability distributions for inter- 
acting and non-interacting cases. (Rb560). 
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FIG. 15. Momentum probability distributions for inter- 
acting and non-interacting cases. (Na530). 



narrow central peak corresponds to atoms localized in the 
lattice sites and the broader background wing to atoms 
which are above the recapture range and do not relocalize 
in the lattice. This resembles the evaporative cooling 
process with narrowed central peak and hot background 
atoms. Cooling here is not dramatic but still present. 
Moreover, the result is in sharp contrast when compared 
to the theoretical and experimental collision studies in 
MOT's where the heating of the trapped atoms due to 
the radiative mechanism is observed pLR but not the 
evaporative-type cooling process. 

The cooling process is observed for all the three atomic 
masses when similar lattice depths are used. The compu- 
tational resources that simulations require, see Sec. VI G, 



do not allow any extensive exploration of parameter space 
but the Csl621 result shows that w ith a deeper lattice 
the situation may change, see Table VII. In shallow lat- 
tices the relative velocity before a collision is small, thus 
enhancing the excitation probability. In deep lattices the 
reduced excitation probability due to large relative ve- 
locity is compensated by the use of more intense lasers. 
The Csl621 result suggests that in deeper lattices one 
may observe heating which is similar to the results from 



MOT studies. But a systematic study of this is out of 
the reach of this paper. 



IX. DISCUSSION AND CONCLUSIONS 

Our results show the basic aspects of one mechanism 
affecting the thermodynamics of the atomic cloud in an 
optical lattice, when the lattice has been prepared with 
near-resonant (detuning a few linewidths), red-detuned 
laser light. In this case the role of inelastic collisions is 
strong, leading to heating and loss of atoms, but this re- 
quires that the interacting atoms are located in the same 
lattice site. This, on the other hand, requires, firstly, 
large atomic densities. Therefore in most lattice stud- 
ies done so far, the role of collisions has been negligible 
due to the low densities, or at least not easy to observe 
(with the exception of collisions producing a clear signal 
such as Penning ionisation of colliding metastable rare 
gas atoms Q). 

Secondly, frequent collisions require clear mobility of 
atoms. This takes place naturally during the Sisyphus 
cooling until the atoms are localized in lattice sites. Thus 
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TABLE VII. Expectation values of kinetic energy per 
atom (< Ek >) for the simulations. The value of p c gives 
the critical momentum which is used in ensemble averaging 
to neglect atoms which have escaped from the lattice. The 
absolute values of the standard deviation are given in paren- 
theses. 
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it is important for suitably dense samples to study the 
role of collisions during the Sisyphus cooling, and our ap- 
proach provides a method which is both dynamical and 
consistent. For the selected parameters our simulations 
show that the Sisyphus cooling process and localization 
of atoms is not prohibited by inelastic processes, i.e., the 
loss and heating of atoms remains small even when the 
average distance between the two atoms is only four lat- 
tice sites. 

Once the localization in lattice sites has been achieved 
as a steady state, the question about the mobility of 
atoms changes to some extent. It should be noted that 
localization does not mean that an atom remains in the 
same site ad infinitum. In the steady state the atoms are 
localized at the sites for most of the time, but also move 
between the sites via tunneling (in the picture where the 
lattice lasers and the excited states are eliminated from 
the effective description). For the selected parameters 
the dipole-dipole interaction does not perturb the lat- 
tice potentials enough to have a significant effect between 
atoms located at different sites (the opposite situation is 
also possible, see Ref. [Q). The tunneling of atoms be- 
tween sites is in the steady state the main process leading 
to inelastic collisions, and as the simulations show (sup- 
ported by the semiclassical estimates) such encounters 
lead mainly to loss of hotter atoms or their selective heat- 
ing. This is because the hotter atoms, naturally, move 
between the sites more frequently than the colder ones, 
creating the selectivity. If the lost atoms are not consid- 
ered, our results show, however, that we can hold on to 
the concept of an existing steady state. 

The collisionally induced velocity-selective loss of hot 
atoms is similar to the evaporative cooling which is 
utilised in magnetic traps to reach ultracold temperatures 
for atoms. It remains to be seen, however, whether the 
resulting cooling effect for the atoms is significant enough 
to be observed in densely filled lattices. In dense samples 
other processes such as reabsorption of photons scattered 
by the atoms are another important source for heating, 
which may well overcome any cooling effect. It should 
be noted, however, that if we ignore the spatial structure 



of the cooling fields, the collisional processes lose their 
velocity-selective nature, and as seen in the simulations 
of Ref. 0, this leads to strong heating of atoms. Here 
the simulations indicate that the lattice structure inhibits 
this heating clearly. 

Our simulations have been very intense computation- 
ally, which makes it very difficult to make the model more 
realistic. Our studies, however, to our opinion, demon- 
strate the basic features to be expected from the colli- 
sions in densely populated near-resonant red-detuned lat- 
tices. There are magnetic-field-assisted cooling schemes 
for blue-detuned lattices for e.g. J = 1 — ► J = 1 systems. 
The blue detuning usually leads to optical shielding, and 
the collisional contribution to inelastic processes is re- 
duced strongly for normal laser cooling intensities as the 
loss channel is expected to be adiabatically closed, see 
Ref. |3^]. As a future prospect it will be interesting to 
study the qualitative differences due to the color of the 
detuning in collisions between atoms in near-resonant lat- 
tices. 
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